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ABSTRACT 

We present numerical radiation-hydrodynamic simulations of cometary H II regions for a number of cham- 
pagne flow and bowshock models. For the champagne flow models we study smooth density distributions with 
both steep and shallow gradients. We also consider cases where the ionizing star has a strong stellar wind, and 
cases in which the star additionally has a proper motion within the ambient density gradient. We find that our 
champagne flow plus stellar wind models have limb-brightened morphologies and kinematics which can see the 
line-of-sight velocities change sign twice between the head and tail of the cometary H II region, with respect to 
the rest frame velocity. Our bowshock models show that pressure gradients across and within the shell are very 
important for the dynamics, and that simple analytic models assuming thin shells in ram pressure balance are 
wholly inadequate for describing the shape and kinematics of these objects at early times in their evolution. The 
dynamics of the gas behind the shock in the neutral material ahead of the ionization front in both champagne 
flow and bowshock type cometary H II regions is also discussed. We present simulated emission-measure maps 
and long-slit spectra of our results. Our numerical models are not tailored to any particular object but com- 
parison with observations from the literature shows that, in particular, the models combining density gradients 
and stellar winds are able to account for both the morphology and general radial velocity behavior of several 
observed cometary H II regions, such as the well-studied object G29.96-0.02. 

Subject headings: H II regions — ISM: kinematics and dynamics — shock waves — stars: formation — stars: 
winds 



1. INTRODUCTION 

Massive stars are born in dense molecular clouds and mod- 
ify their environment by the action of their ionizing photons 
and stellar winds. In a uniform mediu m a spherical H I I 
region, or Stromgren sphere, will form ( Stromgren 1939 ), 



which then expands due to the overpressu re of the pho toion- 
ized gas with respect to its surroundings (Kahn 1954). The 



earliest stages in the lives of high mass stars are generally in- 
visible at optical wavelengths, since the star and its photoion- 
ized region remain embedded in the molecular cloud. These 
clouds are, however, optically thin to radio frequencies and 
infrared radiation and observations at these wavelengths have 
shown that uniform, spherical H II regi ons are the exception. 
The survey at radio frequencies by Wood & Churchwel^ 



for study at all wavelengths show us that the brightest part of 
the nebula corresponds to the ionization front on the surface 
of a molecular cloud, from which material is being photoe- 
vaporated. In the direction away from the cloud the nebula 
is open and ionized gas flows out into the diffuse intercloud 
medium. 

The hydrodynamics of "blister" H II regions was mod- 
eled by Tenorio-Tagle and collaborators in a series of pa 



eled by lenorio-lagle and collaborators in a series or pa- 
pers ( Tenorio-Taglel 1979|; iBodenheimer et alj 19791; Bedijn 



( 198^ developed IRAS color criteria for embedded compact 
and ul tracompact H II regions, which were used in later sur- 



veys ( |Kurtz et aT 1994). These surveys show that around 16% 
of the embedded H II regions have a cometary (i.e., head-tail) 
morphology, while an /RA5-selected methanol maser survey 
shows that 30% of ultr acompact H II regi ons with methanol 



masers are cometaries (Walsh et al. 1998). Around 14% of 



much more evolved, optically visible H II regions, such as the 
Orion Ne bula, also de monstrates cometary shapes at radio fre- 
quencies (Fich 1993 ). In fact, the Orion Nebula gives a clue 
as to the possible origin of the morphology of these objects 
since it is th e archetype o f the so-called "blister" H II regions 
identified by Israel (1978). Its close proximity and availability 
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& Tenorio-Tagle 198 1| ; [Yorke et al.| 1983| ), who renamed them 
"champagne" flows. Later, the effects of including a s tellar 
wind in this scenario were studied by [Comeron ( 1997 ). In 
both of these studies an H II region evolving in a uniform, 
dense cloud breaks out into the diffuse intercloud medium 
once the ionization front reaches the edge of the cloud. Ion- 
ized gas accelerates out of the H II region reaching velocities 
up to ~ 30 km s"'. 

For H II regions still embedded in their natal clouds, the 
classification "compact" is applied to those whose size scales 
are in the range 0. 1-0.3 pc and "ultracompact" to those whose 
sizes are < 0.1 pc. It is thought that these regions, because of 
their small size, represent early stages in the lifetime of the 
star (< 10^ yrs) before the H II region has expanded out of 
the cloud. If a compact H II region is still wholly contained 
within a molecular cloud, a "blister" scenario cannot strictly 
be applied since the ionized region has not reached the edge of 
the cloud. Instead, it could be envisaged that the champagne 
flow occurs due to a density gradient within the cloud. (The 
density gradient could be a relic o f the cloud-collapse pro- 
cess th at produced the massive star) Garcia-Segura & Franco 
( 1996 ) studied the morphologies that arise when ionization 
fronts form in stratified media and spherical density distri- 
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butions, and the instabihties to which they are subject when 
cooHng and stellar winds are included. The hydrodynamics of 
photoevaporate d champagne flows in density gradients was 
investigated by Henney et al. ( 2005 ), who found that long- 
lived quasi-stationary flows are set up in which the ionized 
flow is steady in the frame of the ionization front. 

Observational studies of molecular clouds show that they 
are centrally condensed, with moderate-density envelopes 
(10^-10^ cm"-') surrounding high-density cores (10^ cm"^), 
sugg esting strong density gradients within the clouds (H ofner 
et al. Z000| ). Molecular line studies of dark clouds show 
that the density gradients can be approximated by power-laws 
p{r) cc r"", with e xponents ranging from 1 to about 3 (A r- 
quilla & Goldsmith |1985 ). Near-infrared dust e xtinction stud 



iesof dark clouds show similar results (e.g., Teixeira et al 
2005). Studies of the extended sub-millimeter dust emission 



around ultracompact H II regions also display power-law den- 
sity s tructures, with the exponent v arying between 1.25 and 
2.25 ( Hatchell & van der Tak 2003 ). In these studies, a equal 
to 2 is the most commonly found value for the power-law ex- 
ponent, which corresponds to an isothermal, self-gravitating 
sphere. Analysis of the radio continuum spectra of three 
highly compact H II regions by Franco et al. (200C) also im- 
ply negative power-law density distributions with exponent 
a > 2, possibly even as steep as a = 4 in one case. These 
hyper and ultracompact H II regions are considered to repre- 
sent the earliest stages of H II region evolution, and so provide 
information about the natal cloud density distribution. Once 
the H II region has begun to expand, the density distribution 
of the ionized gas will become shallower as structure within 
the H II region is smoothed out on the order of a sound cross 



i ng tim e. Thus, the density gradients reported by Franco et al 



( pOOCj ) should be treated as lower limits to the initial density 
gradients in the star-forming cores. For the largest values of a 
reported, such a steep radial, power-law density distribution is 
unlikely and it is possible that the density stratification is bet- 
ter represented by a Gaussian or other exponential function. 
Power-law density gradients with exponents 1.6 < a < 2.4 
are a lso reported for giant extragalactic H II regions ( Franco 
et al. |2000[ ), showing that density gradients should be taken 
into consideration at all stages of H II region evolution. 

An alternative explanation for the cometary shape of some 
compact H II regions is that the ionizing star is moving with 
respect to the dense clo i id gas . Mac Low et al. (1991) and 
van Buren & Mac Low ( 1992 ) developed a bowshock mter- 
pretation for cometary compact H II regions, in which the 
ionizing star with its strong stellar wind is moving superson- 
ically with respect to the cloud gas and the ionization front 
becomes trapped in the swept-up shell behind the bow shock 
that forms ahead of the star. Numerical models of supersoni- 
cally moving high mass stars have been carried out by Raga et 



al. ( 1997 ) and L^omeron & Kaper X 1998 ) in the context of run- 
away O stars in the diffuse interstellar medium. Even though 
the parameters are not appropriate for young massive stars 
moving in dense molecular clouds, these numerical models 
show that instabilities complicate the dynamics. Recent nu- 
merical models of H II regions around massive stars orbiting 
in the gravitational potential of structured molecular clouds 
(but without i nclud ing the stellar winds) have been done by 
ranco et al. (2005), who show that a wide variety of unpre- 
dictable morphologies can result. 

Bowshock models are attractive because they easily re- 
produce the limb-brightened morphology seen in radio con- 
tinuum images of cometary compact H II regions, whereas 



champagne flows have no limb brightening. However, it is 
in the kinematics of the ionized gas that these two scenarios 
can really be distinguished. In the bowshock model the high- 
est ionized gas velocities are at the head and are equal to the 
star's velocity. In the champagne flow model the highest ion- 
ized gas velocities are found in the tail and can reach a few 
times the sound speed. Unfortunately, observational kinemat- 
ical studies are difficult to perform given the complete extinc- 
tion that exists in the optical due to the dense environments in 
whi ch these objects are fo und. Radio recombination line stud- 
ies ( IWood & Churchwell||l99l|; lAfflerbach et al.||l994 G aray 



line observations ((^^umsden & Hoare 


1996, 


1999; 


Hoare et al. 


2003; 


Martm-Hernandez et al, 2003 ) 


, and mid-infrared fine- 


structure line observations ( Faffe et al 


. 2003 


; Zhu et al. 


2005) 



have, however, been able to extract some kinematic al infor- 
mation from these regions, which should enable us to distin- 
guish between the various theoretical models that have been 
put forward to explain this phenomenon. It is important to ob- 
tain a good determination of the background molecular cloud 
velocity, which can have a major effect on the interpretation 
of the observations. 

The kinematical observations available suggest that nei- 
ther a bowshock model nor a champagne flow adequately ex- 
plain the velocity variation between head and tail of many 
observed c ometary H I I regio ns, or between the center and 
the wings. 3aume et al. (1994) highlighted the need for more 
realistic models, incorporating stellar winds and more compli- 
cated ambie nt density distributions in th e surrounding molec- 
ular cloud. Lumsden & Hoare ( 1999 ) proposed an empiri- 
cal model for the well-studied object G29.96-0.02 in which a 
champagne flow is modified by the presence of a strong stel- 
lar wind within the H II region which traps the H II region 
in the denser gas into a swept-up shell and diverts the pho- 
toevaporated flow around the hot stellar wind bubble. They 
also required a compo nent of expansio n into the molecular 
cloud of ~ 10 km s"'. Zhu et al. (2005) discuss the need for 
pressure gradients within a paraboloidal shell to explain the 
velocity range for this object. 

In addition to the H II region dynamics, hot, massive stars 
possess strong stellar winds, which will also affect the ve- 
locity structure of these regions. Work by Garcia-Segura & 
Franco ^ 996| ) and preyer et al. ( [2003| ) has shown that the in- 
teraction of the photoionized region with a swept-up stellar 
wind shell produces all sorts of interesting structures, such as 
spokes, shells, clouds and fingers. These arise due to cooling, 
thin shell, Rayleig h-Taylor and oth er instabilities, which have 
bee n analyzed by |Vishniac| ( |l983| ), |Ryu & Vishniac| ( |l988| ) 
and Vishniac ( 1994( ). The shapes of stellar wind bubbles 
forming in H II regio n s with density gradients was studied 
analytically by Dyson (1977). That cometary, compact H II 
regions also contam stellar winds is evident from observations 
of our nearest compact H II region, the Orion Nebula. The 
size scale of the Orion nebula is a few tenths of a parsec in 
diameter Hubble Space Telescope (HST) observations of the 
nebula reveal the presence of bowshocks ahead of the photoe- 
vaporating circumstellar material around the "proplyds" (pro- 
toplanetary disks) in the direction directly towards the main 
ionizing star of the nebula, 0' C Ori. These bowshocks are di- 
rect evidence of an interaction between the photoevaporated 
flows within the nebula and the strong stellar wind from the 
ionizi ng star and enable us to constrain M V,,, of the stellar 
wind dGarcfa-Arredondo et akl^OOlj ^002|). 



In this paper we investigate the formation and evolution of 
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cometary compact H II regions both with and without stellar 
winds in a variety of ambient density distributions and com- 
pare these results to kinematic data from the literature. The 
organization of the paper is as follows: in § | we describe our 
numerical model and test cases, while in § | we describe a 
set of basic reference models comprising a simple champagne 
flow, a bowshock in a uniform medium, and a champagne flow 
with a stellar wind. In § ^ we describe modifications to these 
basic models, including varying the steepness of the density 
gradient and the strength of the stellar wind, and adding stellar 
motion to the champagne flow plus stellar wind models. We 
discuss our results in the context of other models and briefly 
compare our results to observations taken from the literature 
in § Finally, in § ^ we conclude our findings. 

2. NUMERICAL MODEL 
2.1. Method 

We need to calculate both the hydrodynamics and the ra- 
diative transfer within the cometary H II regions. This is a 
non-trivial problem involving a variety of timescales. 

For the r adiative transfer w e adopt a similar procedure to 
that used in Rag aet al.|(|l999|) which uses the method of short 
characteristics ( [Mihalas et al.| [19781 [Kunasz & Aueij|l988[ ) to 
calculate the column density from every source point to each 
point on the numerical grid, and thence the photoionization 
rat e at each p o int. W e have adapted this method, described 



in Raga et al. (1999) for a three-dimensional Cartesian sys- 
tem, to a two-dimensional cylindrically symmetric system. In 
the models presented in this paper we consider only a sin- 
gle source point (i.e., the star) for the ionizing radiation. The 
ionizing flux is calculated from the number of photons with 
energies above the hydrogen ionization threshold of 13.6 eV 
produced by a simple black body of a given effective tem- 
perature and stellar radius. Only one frequency of radiation is 
included in the radiative transfer an d the diffuse field is treated 
by the on-the-spot approximation ( [Henney et al. 2005 ). 

The hydrodynamics is treated with a two-dimensional, 
cylindrically symmetric hybr i d Godunov-van Leer second or 



der m ethod ( |Godunov| |1959| ; [Arthur & Fall^ [1993 



van Leei 



1982). A hybrid method is used in order to avoid the "car- 



buncle" instability seen behind shocks par allel to grid lines 
withpurelyRiemann solver codes (see, e.g., Qniik 1994). Van 
Leer and Riemann solver steps are alternated. The Van Leer 
scheme is no t strict ly upwind according to the definition of 
[Sanders et al. (1998) and puts a small shear term into the mo- 
mentum equation. This shear dissipation is sufficient to over- 
come the carbuncle instability. On the other hand, a Riemann 
solver does a more accurate job in general of representing the 
flow, particularly in transonic regions, and so the combined 
scheme has the benefits of both methods without their respec- 
tive disadvantages. 

The hydrodynamics and radiative transfer are coupled 
through the energy equation, where the source term depends 
on the photoionization heating and the radiative cooling rates. 
The radiative cooling includes collisional excitation of Ly- 
man alpha, collisional excitaton of neutral and ionized metal 
lines (assuming standard ISM abundances), hydrogen recom- 
bination, collisional ionization and bremsstrahlung, with each 
cooling process being characterized by an analytic approxi- 
mation which depends on temperature and ionization fraction. 
The energy equation is solved explicitly using a second order 
Runge-Kutta method. In photoionization equilibrium the gas 
temperature in our models is 10^ K. The time step for the 
calculations is the minimum of the hydrodynamic (Courant 




log(Time) [yrs] 

Fig. 1 . — Expan.sion as a function of time of the ionization front (solid line) 
and preceding shock front (dashed line) for an H II region in a uniform den- 
sity medium. The dotted line represents the analytical expansion law (Spitzer 
1978) for the same parameters. 



condition) and cooling times, evaluated at each time step. The 
present models do not include thermal conduction nor absorp- 
tion of ionizing photons by dust grains. The presence of dust 
within an H II region can have a dramatic effect on the size 
of the Stromgren sphere. However, the distribution of dust 
within the H II region depends strongly on the gas density 
and stellar spectrum and there can even be quite a large vol - 
ume devoid of grains due to sublimation (Arthur et al. 2004). 
Such complications are beyond the scope of this paper. 



2.2. Test cases 

2.2.1. Expansion of an H II region in a uniform density 
medium 

We modeled the expansion of the H II region produced 
by a 30,000 K black body in a uniform medium of density 
6000 cm"^ and initial temperature 300 K. The ionizing pho- 
ton rate of the source is Qho = 2.2 x 10"*** s"' and the theo- 
retical initial Stromgren sphere has a radius of = 0.124 pc 
for an ionized gas temperature of T = 10, 000 K. Our numer- 
ical model also has a Stromgren radius of R„^ - 0.124 pc, for 
an ionized gas temperature of T ^ 10,500 K. Once formed, 
this overpressured sphere of gas begins to expand into the sur- 
roun ding medium. An analytical approximation to the expan- 
sion (Spitzei 1978) is 
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where cu is the sound speed in the ionized gas, Ri is the ion- 
ization front radius and t is time, under a set of simplifying 
assumptions, including (i) that the swept-up neutral shell be- 
tween the shock wave and ionization front is thin such that 
the shock and ionization front move with the same velocity 
(ii) the pressure across the neutral shell is uniform. At late 
times the expansion law is characterized by Ri oc t*^'' . 

In Figure |l] we show the numerical results of a 2D cylindri- 
cally symmetric simulation with grid resolution 512 x 1024 
cells for the shock front and ionization front radii as a function 
of time, together with the analytic result for the same initial 
parameters. At early times the analytic and numerical solu- 
tions differ because the rarefaction wave that propagates back 
through the H II region when the expansion first starts in the 
numerical model reduces the pressure in the ionized gas and 
hence the numerical solution lags behind the analytical one to 
start with. After 10,000 yrs the numerical and analytic results 
agree but then begin to deviate slightly, with the numerical 
solution being ahead of the analytical one at these times. This 
is because the assumptions made for the analytic expansion 
break down. The neutral shell ceases to be thin because the 
shock propagates faster than the ionization front. In particu- 
lar, the density just ahead of the ionization front decreases as 
a fraction of the average density in the neutral shell as time 
goes on, meaning that the ionization front is propagating into 
a successively less dense medium with time, and hence the ex- 
pansion is faster than that predicted by the analytic formula, 
which does not allow for any evolution of conditions in the 
neutral shell. 

2.2.2. Expansion of an H II region in the presence of a stellar 

wind 

The expansion of a stellar wind into a surrounding envi- 
ronment has been studied numerically by man y authors (e.g., 
palleJ|1975|:|Garcia-Segura & Mac Low|[l995|; G ai-cia-Segura 
et al. [l996^ . In one-dimensional studies tlie textbook regions 
of unshocked wind, shoc ked wind, swept-up shoc ked ISM 
are easily distinguished (pyson & Williams] 19971; Weaver 



et al. 1977| ). However, once studies are extended to two di 
mensions, thin shell instabilities in the cooling swept-up ISM 
lead to a more comp lex appearance, which could influ ence the 
H II region outside ([Garcia-Segura & Mac Low| 1995t Garcia- 
Segura et al. |1996|: parcia-Segura & Franco| 1996| ; [Comeron 



1997; Freyer et al [2003). In particular, the knotty appear- 



ance of the unstable thin shell leads to shadowing by clumps 
of the H II region. The resulting pressure variations between 
shadowed and unshadowed regions lead to complicated, low 
velocity (< 2 km s~') kinematics in the ionized gas. 

We used a Chevalier & Cleg j (1985) stellar wind model 
to input the initial stellar wind conditions as a volume- 
distributed source of mass and energy (this method for deal- 
ing with the ste llar wind was also used by |^ozyczka |1985 



1997). This avoids the problem of reverse shocks 



that occur when the stellar wind is input simply as a region of 
high density, high velocity flow near the origin. The mass and 
energy source terms were calculated so as to give the correct 
mass-loss rate and terminal velocity for the stellar wind. The 
volume occupied by stellar wind source terms was chosen to 
be sufficiently small so that the stellar wind always reached 
its terminal velocity and shocked outside of the source region, 
but sufficiently large so as to be a reasonable geometrical rep- 
resentation of a sphere on a cylindrical grid. Of course, this 
thermal wind is not the same as a pressureless stellar wind 
at the center of the inner unshocked wind region, but it does 



have the same ram pressure. The Mach number of the thermal 
wind just before the inner stellar wind shock is ~ 10, and so is 
highly supersonic, meaning that any thermal pressure in this 
region is not dynamically important. 

We performed a number of studies of the evolution of H II 
regions around stellar wind bubbles with the aim of estab- 
lishing the optimum resolution for our calculations. The high 
densities (a; 6000 cm"-') in the ambient gas imply extremely 
short cooling times, requiring a short time step for adequate 
resolution. In addition, a fast stellar wind also requires a 
short time step due to the high temperatures (and hence sound 
speeds) in the shocked wind region. 

On our two-dimensional grids we could not hope to achieve 
the resolution necessary to fully resolve the cooling shocked 
swept-up ISM shell. However, by reducing the cooling rate 
and heating due to photoionization by a given factor we can 
maintain the temperature in the ionized gas while artificially 
enhancing the resolution in the cooling shell (by increasing 
the cooling length). This means that radiative cooling in the 
shell occurs more slowly than it would otherwise have done 
so, but given the densities involved, cooling is still rapid in 
these simulations. We then examined a number of simula- 
tions at different grid resolutions in order to determine a con- 
sistent qualitative behaviour in the cooling region. Note that 
in these simulations there is no need to trigger instabilities 
in the swept -up shell with density fluctuations (e.g., G arcia- 
Segura et al. 1996 ) because the cylindrical symmetry means 
that the stellar wind source volume introduces its own numer- 
ical noise simply by not being perfectly spherical. Further- 
more, we artificially broaden the ionization front (by reduc- 
ing the ionization cross-section), in order to avoid artificial 



nume rical instabilities generated at the grid scale (Williams 



1999) 



In Figure g we show the evolution of the ionized den- 
sity in the region surrounding a (blackbody) star of effec- 
tive temperature Teff = 30, 000 K, stellar wind mass-loss rate 
M - 10 Mo yr ' and terminal velocity y„ - 2000 km s 
where the initial ambient medium has number density no - 
6000 cm"^ and temperature Tq - 300 K. The ionizing pho- 
tons and the stellar wind are started simultaneously. The top 
four panels show the ionized density at times 200, 500, 1000 
and 1500 years since the start of the calculation with grid res- 
olution 400 X 400 cells. The grid shown represents a spatial 
size of 0.15 X 0.15 parsecs. The H II region forms essentially 
immediately, on a recombination timescale, with a radius of 
0.13 parsecs. The stellar wind bubble takes longer to form and 
will grow until balance is achieved between the ram pressure 
of the stellar wind and the thermal pressure of the surrounding 
ISM. 

By the time of the first panel (200 years), instabilities have 
begun to form in the swept up shell of ionized gas, at 0.03 par- 
secs radius. The rapid cooling in this shell means that it is 
t hin and subjec t to the thin shell and Vishniac instabilities 
(|Vishniac||l98"3|; |Ryu & Vishniac| |l987|, |1988| ; [Vishniac & Ryu 



1989; Vishniac 



1994 ), which have been shown to occur in 2D 



num erical stellar wmd bubble models (Garcia-Segura & Ma c 



19951 ; parcia-Segura & Franco||1996| ; [Freyer et al.|[2003| ) 



Low 

The ionized densities in the shell reach lO"" cm"-', i.e., com- 
pression factors much higher than the factor 4 expected for 
an adiabatic shock, confirming that this is a region of strong 
cooling. The increased opacity due to the swept-up shell leads 
to recombination of the ionized region beyond. The inhomo- 
geneities in the shell represent opacity variations, which lead 
to variations in the recombination rate that can be detected as 
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Fig. 2. — Evolution of logio(density) in an H II region with a stellar wind. Grayscale for the ionized gas covers the range 70 < n, < 70000 cm ^, while for the 
neutral gas the grayscale range is 200 < n„ < 200000 cm"', (a) Ionized gas density after 200 years (full range < n, < 11700 cm"'), showing the formation of 
ripples in the swept-up thin shell due to cooling instabilities and corrugations in the ionization front due to shadowing, (b) Ionized gas density after 500 years 
(0 < itj < 92400 cm"') showing the growth of the shadowing instability, (c) Ionized gas density after 1000 yeai's (0 < n; < 61 100 cm"'), (d) Neutral gas density 
after 1500 years (0 < n„ < 210000 cm"') showing the formation of dense neutral clumps at the base of the ionized spokes, (e) Same as (c) but at half the grid 
resolution (0 < «, < 45200 cm"'), (f) Same as (e) but with a lai'ger cooHng coefficient (0 < n, < 63000 cm"'). 



ripples in the outer rim. 

In the second panel (500 years), the instability has saturated 
and the knotty thin shell at 0.045 parsecs is being pushed out- 
wards by the hot, shocked stellar wind bubble. The degree of 
ionization in the H II region beyond the swept-up shell is very 
small as this region has now essentially completely recom- 
bined. The knotty nature of the thin shell causes variations 
in the optical de pth and th i s lead s to the shadowing instabil- 
ity discussed by [Williams| ( 1999 ). Spokes of ionized gas are 



characteristic of this instability and represent the propagation 
of R-type ionization fronts in collimated beams through gaps 
(regions of underdensity) in the knotty shell. The bases of 
the beams are wider than the tips because the ionized gas at 
the base is pressurized and the gas can expand laterally. The 
freely flowing and shocked stellar wind regions interior to the 
swept-up shell are of such low density that they do not con- 
tribute to the map of ionized density despite being fully ion- 
ized regions. 

In the third (1000 years) and fourth (1500 years) panels 
the knotty shell has expanded outwards pushed by the hot 
shocked stellar wind, reaching a mean radius of 0.08 parsecs 
in panel four. The dense clumps of neutral gas that form at the 
base of the ionized spokes are the result of the shadowing in- 



stability ( [William s||1999| ). Their density continues to increase 
as time goes on, unlike the ionized gas density in the swept- 
up thin shell, which decreases since the pressure in this region 
decreases as the expansion proceeds (the temperature of this 
gas remains constant at 10"* K). The main ionization front is 
trapped in the swept-up stellar wind shell at these times. As 
the expansion of the stellar wind bubble proceeds, the move- 
ment of the knots, plus the contributions to the opacity of the 
material photoevaporating from the neutral clumps, lead to 
changes in the pattern of spokes. 



The shadowing instability discussed by Williams ( 1999 ) re 



suits from the effect on an R-type ionization front of inhomo- 
geneities in the upstream density. First, corrugations in the 
ionization front are formed, which become amplified by the 
instability, leading to the formation of characteristic spokes of 
ionized gas with clumps of dense, neutral material at the base. 
In the simulations presented in Figure^, the shadowing effect 
is produced by clumps formed by the thin-shell instability in 
the swept-up stellar wind shell, within the H II region itself. 
Although corrugation of the outer ionization front is seen at 
early times (Fig. ^), this is not important for the growth of 
the instability, since the whole region outside of the thin shell 
is recombining by this point. The shadowing instability takes 
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TABLE 1 
Model parameters. 



Model 
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H 


M 




V, 




n 




pc 




km s"' 


km s"' 


A 


no expiz/H) 


8000 


0.05 








B 


no 


6000 




10-* 


2000 


20 


C 


no exp{z/H) 


8000 


0.05 


10-'^ 


2000 




D 


no exp{z/H) 


8000 


0.05 


10-^ 


2000 




E 


noexpizlH) 


8000 


0.2 


10-'' 


2000 




F 


«o[l + («/0.01)2]"' 


9.8 X lO*" 




lO-*" 


2000 




G 


no expizlH) 


8000 


0.05 


10-*^ 


2000 


5 


H 


no exp(zlH) 


8000 


0.05 


10-'^ 


2000 


10 


I 


no exp{z/H) 


8000 


0.2 


lo-* 


2000 


10 



hold within the thin shell itself, where the main ionization 
front has become trapped. Because the whole structure is ex- 
panding due to the high pressure of the shocked stellar wind, 
the pattern of ionized spokes in our simulations changes with 
time. 

The fifth panel of Figure || shows the effects of reducing the 
grid resolution. Reducing the resolution by a factor 2 in each 
dimension gives qualitatively similar results to the standard 
case. 

In the sixth panel we have used the full value of the cooling 
rate and heating rate with the lower grid resolution. Increas- 
ing the cooling and heating rates means that cooling occurs 
sooner and more quickly, leading to a thinner, denser swept- 
up shell. The recombining region has receded to 0.12 parsecs 
by 1000 years of evolution while the swept-up wind shell is 
located at 0.05 parsecs from the central star, as compared to 
0.06 parsecs in the standard case. However, for this grid res- 
olution the shell is certainly not resolved and would require 
much higher resolution. In these calculations it is desirable to 
resolve the shell as well as possible because of the effect this 
has on the ionized region beyond. The device of decreasing 
the cooling and heating rates permits us to better resolve the 
swept-up shell without having to use an unpractical number 
of grid cells. 

3. BASIC MODELS 

Our models are motivated by a desire to explain features 
of the spectra of cometary ultracompact H II regions such as 
G29.96-0.02, whose morphology and kinematics cannot sat- 
isfactorily be explained by either a simple champagne flow 
or a simple bowshock model. In particular, we are interested 
in explaining the turnover from blueshifted to redshifted to 
blueshifted velocity seen going from the head to the tail along 
the symmetry axis in observat ions of this object in hydrogen 
recombination lines (see , e.g., Lumsden & Hoare 1996, 1999 



[Vlartm-Hernandez et alj|2003 ) and fine-structure lines of the 
Ne^ ion OZhu et al.| |2005[ r 

Instead of tailoring a model to fit a particular observation, in 
this paper we wish to identify features in the simulated spectra 
that correspond to different configurations of density distribu- 
tion, stellar wind parameters and stellar motion. To this end, 
we consider only one set of parameters for the ionizing star, 
namely a 30,000 K black body with an ionizing photon rate 
of 2.2 X 10 ^^ s~'. This co r responds approxima tely to spectral 



type 09V ( |Panagia||1973| ; [Martins et al.|[2005D . The ambient 
medium is taken to be at rest and neutral, with a temperature 
of 300 K. For the majority of models, the grid resolution is 



750 X 1000 cells, representing a total grid size of 0.6 x 0.8 par- 
secs. 

To begin with we discuss three "canonical models". Mod- 
els A, B, and C of Table |lj which correspond to a simple 
champagne flow, a stellar bowshock in a uniform medium, 
and a champagne flow with a stellar wind, respectively. The 
model parameters are summarized in Table ll|. Model A is fol- 
lowed for 20,000 yrs, while Models B and C are evolved for 
40,000 yrs. 

3.1. Model Parameters 

3.1.1. Champagne flow: Model A 

We set up a photoevaporated flow in an exponential density 
distribution, where the density at the position of the source is 
no - 8000 cm"^ and the scale height is H - 0.05 pc. 

3.1.2. Bow Shock: Model B 

The source is assumed to move with a velocity of 20 km s"' 
in a uniform medium of density «o - 6000 cm"-'. This veloc- 
ity is mo tivated bv previous empirical rnodels (e.g., A ffler- 
bach et al. |l994|; [Lumsden & Hoai-e||l996|; Martin-Hernandez 
et al. [2003t [Zhu et al.[ [2005[ ) for the ultracompact H II re- 
gion G29.96-0.02. This star has a strong stellar wind with 
mass loss rate M - 10"'' Mq yr"' and terminal velocity 
Vw - 2000 km s"', typical of an early type star. After the 
20,000 yrs of evolution followed by the simulation, the star 
has traveled a distance of 0.4 pc. The resolution of this model 
is 5 12 X 1024 cells, representing a spatial size of 0.4 x 0.8 par- 
secs. 

3.1.3. Champagne flow with stellar wind: Model C 

We investigate the eff'ect of including a stellar wind in the 
champagne flow described in Model A. The stellar wind has a 
mass loss rate M - IO^^Mq yr"' and terminal velocity - 
2000 km s"', i.e., the same wind parameters as used in the 
standard bowshock model. 

3.2. Basic model results 

3.2.1. Velocity field 

Figures H ^ and H show the logarithm of the ionized density 
(greyscale) and velocity field in the r- z plane for each of the 
canonical models A (champagne flow), B (bowshock), and C 
(champagne flow plus stellar wind) described in Table |l] after 
10,000 yrs, 20,000 yrs and (for Models B and C) 40,000 yrs 
of evolution. The velocity field is represented by a split scale: 
the black arrows show the highest velocity gas y > 30 km s"' 
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Fig. 3. — Logarithm of ionized density for a cut thi'ough tlie r - z plane of 
the numerical simulation of Model A (simple champagne flow) after elapsed 
times of 10,000 yrs (top panel) and 20,000 yrs (bottom panel). The arrows 
show the velocity field, with the scale arrow representing 30 km s"' in these 
models. The position of the source is marked by a black triangle. 



(generally representing stellar wind gas), while the white ar- 
rows are scaled to the lower velocity gas u < 30 km s"' 
(present in the photoevaporated flows and swept up shells). 
The grayscale shows the densest ionized gas in black with a 
dynamic range of 2 orders of magnitude. Neutral gas, or very 
low density ionized gas (i.e., stellar wind) appears white. 

From these figures we see that in the case of the pure cham- 
pagne flow. Model A, the ionization front is smooth, whereas 
that of Model C, which has a stellar wind, has a rippled ap- 
pearance due to the formation of instabilities in the swept-up 
wind shell, such as those presented in Figure ^ and discussed 
above. In Model B, the ionization front becomes trapped in 
the swept-up wind shell, which forms a bowshock ahead of 
the moving star. In Model C the ionization front is trapped 
ahead of the star in the denser gas but escapes at the sides 
where the density is lower and forms a champagne flow here. 
As regards temporal evolution, we see from Figure || that 
the champagne flow model (Model A) does not appreciably 
change in appearance between the two times depicted . This 



quasi-stationary phase was described by Henney et al. (2005) 



and the ionization front advances very slowly (~ 2 km s"') 
up the density gradient. At the earlier time, the champagne 
flow has not completely cleared the low-density material from 
the grid: the rarefaction wave moving back from the initial 
Stromgren surface through the ionized gas when the H II re- 
gion begins to expand hydrodynamically reaches the symme- 
try axis earlier where the ionization front shape is narrower 
This shows up as a "bump" in the velocity distribution of the 
upper panel of Figure |[ By the time 20,000 yrs have passed, 
the rarefaction wave has reached the symmetry axis at all z 
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Fig. 4. — Same as Figure ^but for Model B (simple 20 km s"' bowshock 
model) and showing evolution after 10,000 yrs (top panel), 20,000 yrs (mid- 
dle panel) and 40,000 yrs (bottom panel). The black arrows represent veloci- 
ties of 30 < I) < 2000 km s" ' in this figure, while the white arrows represent 
velocities of « < 30 km s" ' . 



heights and the champagne flow fills the whole grid. We note 
that for this m odel w e do not see the shadowing instability 
discussed in § 2.2.2 because the initial conditions are totally 
smooth and there is no stellar wind shell to form clumps. 

In Figure ^ we can see how the presence of a strong stellar 
wind affects a champagne flow. The initial ionization front 
distance ahead of the star for Models A and C is the same, 
reflecting the fact that the ionized region is set up within 
one recombination time of the star switching on. The ion- 
ization front distance in Model C then retreats, which corre- 
sponds to the formation of the dense swept-up stellar wind 
shell that traps the ionization front. At these very early times 
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Fig. 5. — Same as Figure ^ but for Model C (champagne flow plus strong 
stellar wind). 



(f < 4000 yrs), before the density gradient has become im- 
portant, the shadowing instabihty produces spokes of ionized 
gas protruding beyond the swept-up shell. However, once the 
bubble diameter is greater than about one density scale height, 
the strong advection around the shell due to the density gra- 
dient means that any clumps that could cause shadowing are 
moving too fast for the instability to take effect. The initial 
shadowing does open up "holes" in the thin shell, through 
which the hot, shocked stellar wind can flow, generating a tur- 
bulent interface. As the hot, pressurized, energy-driven stellar 
wind bubble expands supersonically, the ionization front ad- 
vances up the density gradient and after 5000 yrs has caught 
up with the simple champagne flow ionization front. Even- 
tually, in the direction of decreasing density the stellar wind 
bubble blows out, and there is a transition to a momentum- 
driven flow. Meanwhile, a champagne flow starts to establish 
itself within the swept-up shell. Initially, the champagne flow 
is blocked from expanding down the density gradient because 
the stellar wind is strong enough that the contact discontinu- 
ity between the shocked wind and the swept-up material is 




2x10^ 



3x10' 



Time [yrs] 



Fig. 6. — Distance from star of ionization front propagating up the density 
gradient against time for a champagne flow (thick solid line) and a champagne 
flow plus stellar wind (thick dotted line). Also plotted are the positions of 
the respective shock fronts in the neutral material ahead of the ionization 
front: thin solid line — champagne flow, thin dotted line — champagne flow 
plus stellar wind. 



pushed past the position where the champagn e flow would 
go th rough a sonic point (see, e.g., case A of Henney et al 



EOOS) and so the champagne flow stays subsonic on the axis 
The champagne flow first establishes itself off-axis, where the 
advancing ionization front opens up and distances itself from 
the contact discontinuity, allowing the flow to become super- 
sonic. The opening up of the ionization front appears to be 
driven by Kelvin-Helmholtz instabilities, which form at the 
head and travel down the structure, broadening as they reach 
regions of lower density and pressure. This leads to varia- 
tions in the opacity and allows the ionization front to expand. 
Even though the external shock in Model C is the stellar wind 
shock, rather than the ionization front shock of Model A, the 
shock velocity is so small (see Figure ^ that the density en- 
hancement is negligible and so the ionization front in Mod- 
els A and C follows the same trajectory. The majority of the 
opacity in this part of the H II region comes from just behind 
the ionization front, since this is where the densest gas is to 
be found, hence even though the stellar wind has evacuated a 
cavity behind the ionization front, the conditions in the pho- 
toionized gas behind the ionization front are almost the same 
as in the pure champagne flow case and so the expansion is 
the same. 

In Model B, at very early times, the shadowing instability 
is important (not shown here) and is triggered by a coolin g 
instability in the swept-up thin shell (as discussed in § 2.2.2 ), 
which produces shadowing clumps. At the base of the spokes 
of ionized gas, the gas pressure widens the "gap" in the shell 
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and allows the shocked stellar wind to flow through. This can 
be seen in the top panel Figure ^ as "ears" on the bowshock 
structure. Although strong advection around the bowshock 
after f > 10, 000 yrs leads to the disappearance of the charac- 
teristic spokesof ionized gas, the initial instability is the seed 
for a very turbulent interface between the hot, shocked stellar 
wind and swept-up material. This turbulent interface leads 
to the shedding of large eddies from the head of the bow- 
shock and affects the thickness of the trapped ionized region. 
The eddies contain large amounts of dense, ionized material. 
At later times, Kelvin-Helmholtz instabilities at t he interface 
region seem to be the domin ant instability (cf. |^aga et aT 



1997; I^omerdn & Kapei 199a). The time sequence of images 
shown in Figure ^ shows that the shape of the bowshock varies 
considerably. At 10,000 yrs the shell is very thin, but with a 
large bowshock opening angle, at 20,000 yrs the shell is thin 
at the head but thicker at the sides, with a narrow bowshock 
opening angle, while at 40,000 yrs, the shell is quite thick and 
the opening angle of the bowshock is also quite large. The 
dense "blobs" on the axis in the 20,000 yrs figure are a result 



of the imposed cylindrical sy mmetry focussing the turbulent 
flow onto the symmetry axis. Raga et al. (1997) state that, for 
late times of their models of high-velocity runaway O stars 
in low-density media, the bowshock flow reaches a "statis- 
tically stationary" regime. Our models have yet to achieve 
this regime but it could be expected that they will, since by 
40,000 yrs all the transient features due to the initial insta- 
bilities have been advected a long way downstream. The ed- 
dies that are shed periodically from the apex of the bowshock 
cause pressure variations that alternately pinch then open up 
the region of shocked stellar wind. 

The velocity fields in these models are quite complex. The 
photoevaporated flows in the champagne models (A and C) 
start at lo w velocities at the ioni zation front (v ~ 5 km s"', see 



case A of Henney et al. (2005)) and are initially accelerated 
rapidly to the sound speed (y ~ 10 km s '), and then more 
gradually due to the density gradient, reaching velocities of 
around v ^ 30 km s"' before leaving the numerical grid (this 
gas would be below the detection limit due to its low density). 
In Model C, with a stellar wind, we also see thin shell insta- 
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Fig. 8. — Simulated emission line moments for Model A (simple champagne flow), inclination angle 45° — line intensity, mean velocity and RMS velocity width 
(multiplied by scale factor 2 yj2 log 2) as a function of position along the slit. Left: sHts parallel to symmetry axis; Right: sHts perpendicular to symmetry axis. 
The offsets of the slits from the star are indicated in the upper left corner. 
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Fig. 9. — Same as Fig. |^ but for Model B (20 km s ' bowshock). 
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bilities, which form first in the densest part at the head of the 
cometary-shaped thin shell (because of faster cooling here) 
and later round the sides. The resulting knots are then de- 
formed both by advection within the swept-up shell and also 
by Kelvin-Helmholtz instabilities due to the shear between the 
shell flow and the photoevaporated flow. Large knots of dense 
ionized material formed at the head of the cometary flow can 
thus be advected off the grid during the simulation timescale. 

In the bowshock case (Model B), velocities at the head are 
roughly equal to the stellar velocity. Towards the tail of the 
cometary H II region, velocities in the dense, ionized gas are 
low and randomly oriented due to eddies formed at the shear 
interface. Although the highest velocities are to be found in 
the stellar wind itself, this gas has such low density that it is 
not detectable. 

3.2.2. Simulated emission measure maps 

Figure ^ shows the simulated emission measure smoothed 
by a Gaussian beam of 1.0". The numerical grid was pro- 
jected onto a pixel grid whose longest dimension of 80 pixels 
represents a spatial size of 2R cos / -i- Z sin /, where R and Z 
are the maximum extents of the radial, r, and longitudinal, z, 
directions of the numerical model (see e.g.. Figure 0), and / is 
the inclination angle to the line of sight. In Figure^ we have 
arbitrarily assigned an angular size of 40" to be equivalent to 
80 pixels. This would place all Models A and C at distances of 
approximately 7 kpc, depending slightly on inclination angle, 
while Model B would be at approximately 6 kpc. In Figure^ 
we plot only a 24" x 24" subset, with the projected stellar 
position assigned to (0, 0) in each case. For reasons of space, 
we show results only for inclination angles / = 30°, 45° and 
60°. Emission measure maps for other inclination angles are 
available from the authors on request. 

The emission measure maps show a range of morpholo- 
gies. Model A (champagne flow) gives a clump or fan-shaped 
emission, as expected, while the addition of a stellar wind 
(Model C) flattens this emission and gives a more arc-like ap- 
pearance. The bowshock model (Model B) results in an arc of 
limb-brightened emission, as expected. 

3.2.3. Simulated spectra 

In Figures &-[lQ we plot the line intensity, mean velocity 
and velocity dispersion for slits parallel and perpendicular to 
the symmetry axis for the models shown in Figure ^ for incli- 
nation angle 45°. Results for other angles are available from 
the authors on request. These quantities are evaluated using 
the emission line moments, where the velocity dispersion is 
the RMS width scaled by a factor 2 -y/2 log 2, which would 
give the FWHM in the case of a Gaussian profile. In each 
case the slit width is one arcsec and we have included thermal 
broadening using the temperatures which result from the nu- 
merical simulation. In the gas in photionization equilibrium 
this temperature is 10"* K but the partially ionized gas (either 
recently photoionized or in the process of recombining) will 
be at a lower temperature, and gas in the ionized, shocked 
stellar wind will be at a much higher temperature. The flux 
intensity for each slit is relative to the maximum intensity of 
the whole object. The slit offsets indicated are relative to the 
symmetry axis (for the parallel slits) and stellar position (per- 
pendicular slits). 

The peak flux intensity occurs on the symmetry axis ahead 
of the star. For the slits parallel to the symmetry axis, the flux 
intensity decreases steeply along the slit for all models. For 



the slits perpendicular to the symmetry axis, there are differ- 
ences between Models A, B, and C. For Model A (champagne 
flow), the peak flux in each slit always occurs on the symmetry 
axis and falls off to the sides. For Models B (bowshock) and 
C (champagne flow plus wind), the flux intensity has a flat- 
topped shape as the slits move down the tail of the cometary 
H II region, with peaks at the sides due to the emitting mass 
being concentrated in a shell around the low-density stellar 
wind bubble. 

The velocity structures of Models A, B, and C are very dif- 
ferent. In Figure ^we see that the flow in Model A leaves the 
head of the cometary H II region with blueshifted velocities. 
It then accelerates and reaches maximum (blueshifted) speed 
at the projected position of the star Thereafter, the velocities 
become more redshifted because of the contribution from the 
higher density, lower velocity, redshifted gas on the other side 
of the object in projection. 

In Figure ^ (Model B) the velocities at the head are red- 
shifted and correspond to the bowshock region being pushed 
forward ahead of the moving star. The noisiness of this spec- 
trum is due to the eddies that form in the shear flow inter- 
face between the shocked stellar wind and the swept-up am- 
bient medium. The drop in velocity just ahead of the pro- 
jected position of the star in the slits parallel to the symmetry 
axis is because the line-of-sight cuts through a large region 
of low-density, low-random-velocity gas as well as through 
the more ordered, high blue-shifted velocity, dense material 
in the swept-up shell in the bow region. The velocities fall 
off steeply towards the tail region because the influence of the 
stellar motion on the gas is very small by this point. 

Figure 1^ (Model C) shows a combination of blue-shifted 
and red-shifted material. As for the simple champagne flow 
(Model A), the photoevaporated flow leaving the ionization 
front is blue-shifted. However, the stellar wind confines the 
champagne flow into a shell, causing it to flow around the 
stellar wind bubble. The line of sight through the shell at 

offset 7" cuts through a thick, dense region of redshifted 

velocities, which weight the spectrum and lead it to turn over 

The velocity widths of Model A increase from head to tail 
of the champagne flow but are uniform across it for a given 
slit, and do not differ much from the value of the thermal 
width. In Model B, the velocity widths are noisy for slits 
parallel to the symmetry axis, due to the eddies formed in 
the shear flow. The widths of Model C increase towards the 
wings of the cometary H II region, due to the ionized gas be- 
ing channeled into a shell and accelerated around the stellar 
wind zone. 

Although we do not show them here, the spectra for Mod- 
els B and C at 40,000 yrs are very similar to those at 
20,000 yrs, particularly as regards the velocity range and pat- 
tern. There are differences in the line widths, with the widths 
after 40,000 yrs being generally smoother, reflecting the fact 
that eddies formed earlier in the simulation at the head of the 
object are advected down the structure where they expand and 
so their influence diminishes. 

4. MODIFIED MODELS 

Taking Models A, B, and C as our reference models, we 
vary parameters such as the stellar wind strength and the den- 
sity scale height to examine their effect on the morphology 
and kinematics of the resultant cometary H II regions. We 
also study the effect of stellar motion up a density gradient, 
considering stellar speeds of 5 and 10 km s The parameter 
space to explore is huge, and we restrict ourselves to a few 
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representative cases. 

4.1. Model parameters 

4.1.1. Champagne flow plus stellar wind: Models D, E, and F 

These models are similar to Model C, except that in Model 
D the stellar wind mass loss rate is an order of magnitude 
smaller, M - IO^^Mq yr"', in Model E the scale height of the 
density distribution is much larger, H - 0.2 parsec, and hence 
the density gradient is shallower In Model F the density dis- 
tribution is spherical, rather than plane parallel, and is given 
by 

1 + 



no 



(2) 



where no = 9.8 x 10'' cm"^ is the number density at the center 
of the core of this distribution, where the core is assumed to 
have a radius r^ - 0.01 pc, and - r^ + is the spherical 
radius. The exponent o- = 1 in this model, giving n oc R-^ at 
large radii. The star is offset from the center of the distribution 
such that n - 8000 cm"-' at its position. 

4.1.2. Champagne flow plus stellar wind plus stellar motion: 
Models G, H, and I 

These models combine a density gradient, stellar wind and 
stellar motion. Models G and H have an exponential den- 
sity distribution with scale height H - 0.05 parsec. In 
Model G the star's velocity is 5 km s"' and the star starts 
from 0.05 pc below the position where the number density is 
n = 8000 cm"^, while for Model H the velocity is 10 km s"' 
and the star starts from 0. 1 pc below this reference position. 
Finally, Model I also has an exponential density distribution 



but with a larger scale height, H - 0.2 pc and the star's veloc- 
ity is y* - 10 km s"', where the star again starts its motion a 
distance 0. 1 pc below the position where the number density 
is n = 8000 cm-\ 

4.2. Modifled model results 

4.2.1. Velocity fleld 

Figure [ill shows the ionized density (greyscale) and veloc- 
ity field in the r - z plane for Models D, E and F of Table |l] 
after 20,000 yrs of evolution. Model D looks very similar 
to a champagne flow, with the stellar wind only affecting a 
narrow region deep within the ionized zone. Models E and 
F are both cases of a champagne flow in a shallow density 
distribution. In these cases, the stellar wind can reach pres- 
sure balance with the external medium within a scale height 
and the inner wind shock is almost spherical. The shocked 
stellar wind region is also approximately spherical, except for 
in the tail region where it has broken out into the lower den- 
sity surroundings. At the head of the cometary H II region 
there is a thick shell of ionized gas flowing off the ionization 
front, which is then diverted around the outside of the stel- 
lar wind bubble. Because the ambient density does not fall 
off" very quickly, the density in the ionized gas shell is high 
and fairly uniform, and the ionization front becomes trapped 
in the shell, unlike in Model C, where the rapid fall-off" in 
density means that the ionization front can break through the 
swept-up shell in the tail. In Models E and F the shallow den- 
sity gradient allows the shadowing instability to form broad 
spokes in the lowest density gas, where the pressure at the 
base of the spokes opens gaps in the swept-up shell through 
which the shocked stellar wind can flow. This can be seen 
clearly in the corresponding panels of Figure |ll[ 
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Fig. 1 1 . — Same as Fig. ^ but for models (from top to bottom) (a) Model D: 
champagne flow plus stellar wind, H = 0.05 pc, M = IO^'Mq yr"' (b) Model 
E: champagne flow plus stellar wind, H = 0.2 pc, M = IO^^Mq yr"' (c) 
Model F: champagne flow, spherical density distribution. All models shown 
after 20,000 years of evolution. 



The velocities in Model D are very similar to the simple 
champagne flow model. In Models E and F the velocities in 
the ionized shell are fairly low and there are some turbulent 
regions due to instabilities which are slowly being advected 
down the shallow density gradient. 

Figures [l^ and [l^ show the ionized density and velocity 
fields for Models G, H, and 1. Model 1 has the shallowest den- 
sity gradient and this is reflected in its more spherical, thicker, 
ionized shell. There is a lower emission measure region out- 
side of the main shell. Models G and H have the same steep 
density gradient. Both of these models show a champagne 
flow as well as a stellar wind swept-up shell. The stellar ve- 
locity in Model G is half that of Model H and this results in the 
photoionized shell ahead of the star in Model G being thicker 
than that in Model H. In Figure |l2|, the star in Model H has 
reached a denser part of the ambient density distribution than 
the star in Model G but comparing the two models when the 
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Fig. 1 2. — Same as Figure ^ but for models with density gradient plus 
stellar wind plus stellar motion, (a) Top — Model G: V, = 5 km s"', H = 
0.05 parsec. (b) Bottom— Model H: V, = 10 km s"', // = 0.05 parsec. Both 
models shown after 20,000 yrs of evolution. 
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Fig. 13. — Same as Figure ^but for models with density gradient plus stellar 
wind plus stellar motion (from top to bottom) (a) Model I: V» = 10 km s"', 
H = 0.2 parsec after 10,000 yrs (b) Model 1 after 20,000 yrs of evolution. 
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star is at the same position still shows that the photoionized 
shell ahead of the star in Model H is thinner Also, the hot, 
shocked stellar wind in Model G occupies a greater fraction of 
the volume than in Model H. In fact, in the tail of Model H, the 
stellar wind does not go through a global decelerating shock 
as in Model G, but instead appears to have an oblique shock 
surface that neither heats nor decelerates the wind very much. 
In Model 1, the stellar wind suffers a global decelerating shock 
but then in the downstream direction is focussed by the ion- 
ized shell and reaccelerated by a nozzle effect. 

The velocity fields of Models G, H, and 1 show that large 
accelerations can occur around the sides of the flow in the 
ionized gas, with velocities reaching 30 km s"' or more. This 
is an effect of including a density gradient, since in a simple 
bowshock model there is no acceleration around the sides and 
velocities fall off quickly towards the tail. Models G and H 
show the largest accelerations towards the tails of the cham- 
pagne flows, and this is entirely consistent with them having 
the steepest density gradient. In Model I, the highest veloc- 
ities occur in the swept up shell as it expands into the low 
density region. Ahead of the star, the velocities at the ioniza- 
tion front in these three models are low (i.e., lower than the 
stellar velocities) and in the direction of stellar motion, unlike 
the stationary star cases. 

In these models, the shadowing instability only occurs at 
the earliest times and is most noticeable in Models G (due 
to the low stellar velocity) and 1 (due to the shallow den- 
sity gradient). In both cases the instability soon disappears 
{t < 3000 yrs) because of advection of the shadowing clumps 
around the shell. At later times the dominant instabilities are 
Rayleigh-Taylor and Kelvin-Helmholtz type due to the shear 
between the champagne flow and the shocked stellar wind 
flow. 

4.2.2. Emission measure maps 

In Figure |l^ we show the emission measure maps at pro- 
jection angles of / = 30°, 45° and 60° for models D, E, and F. 
(Emission measure maps for other inclination angles are avail- 
able from the authors.) Model D closely resembles the emis- 
sion measure one expects from a pure champagne flow, since 
the stellar wind is not strong enough (M = lO^^M© yr ') to 
confine the photoevaporation flow from the ionization front. 
Subtle evidence of the wind interaction can be seen in the 
contours. There is a lack of emission on the axis and a faint 
brightening corresponding to the shell formed by the interac- 
tion between the wind and the champagne flow. 

Models E and F are examples of champagne flows with stel- 
lar winds in shallow density gradients. In these cases, the 
stellar wind is strong enough (M = IO^^Mq yr"') to con- 
fine the photoevaporated flow to a thick shell and the resul- 
tant emission measure maps show an arclike rather than fan- 
like morphology. The shallow density gradient in these mod- 
els means that the emission measure does not fall off very 
quickly towards the tail of th e H II region and so they ap- 
pear more shell-like (see, e.g., De Pree et al. 2005 ). Models 
G, H, and I all show limb-brightened emission measure maps 
(Fig. 15), with the limb-brightening being most enhanced for 
the largest angle to the line of sight (/ = 60°). Model I, having 
the shallowest density gradient, has a "filled shell" type mor- 
phology, whereas Models G and H are arclike. In Model H 
the arc of most intense emission is situated much closer to the 
star than in the other models, principally because the star has 
reached denser gas than in Models G and I. The ionized shell 
in Model H (see Figure n2h is much thinner than in Models G 



and I and this is reflected in the steep drop-off in intensity both 
ahead and behind the emission peak. 

4.2.3. Simulated spectra 

The spectra of Model D shown in Figure [l^ do not differ too 
much from the standard champagne flow spectra (Model A). 
In the spectrum of Model E (Fig. [Tt]) the flux does not fall off 
so fast with distance from the headof the cometary H II region 
since the density gradient is shallower. Knots in the swept- 
up wind shell cause features in the velocity plot, but other- 
wise line-of-sight velocities in the ionized gas are very low (< 
5 km s '). The density gradient has little effect on the dynam- 
ics because it is so shallow and the swept-up shell of ambient 
material is neither expanding very rapidly nor has important 
transverse motions. The spectrum of Model F (Fig. |l^) is 
very similar, although for offsets between -10" and -5" the 
flux, velocity and velocity widths are dominated by the bro- 
ken shell which can be seen in Figure This broken shell 
is a transient feature which will eventually be advected down- 
stream. All the Models G, H, and I show blueshifted ve- 
locities at the head of the cometary H II region changing to 
redshifted velocities once the offset has passed the position of 
the star, as can be seen in Figures p^-|2l|. Even though the 
ionized gas velocities ahead of the stars in these models are 
positive (i.e., in the direction of stellar motion), the projection 
angle gives more weight to the gas in the blue-shifted flow in 
the ionized shell and champagne flow at the nearside of the 
cometary H II region. The "noise" seen in these spectra is 
due to instabilities (eddies) moving through the ionized layer 
Model I, having the shallowest density gradients, shows the 
smallest range of velocities (~ 10 km s"'), while models G 
and H show a range of ~ 15 km s ' even though the velocity 
of the star is only 5 km s"' and 10 km s"' respectively. This 
is mainly due to the champagne flow formed in the steep den- 
sity gradient and is very similar to what is seen in the case of 
Model C (Fig. |l^. 

The velocity widths show varied behaviour. Model G 
achieves widths of ~ 40 km s"' in the wings of the tail. On 
the other hand, the highest widths ~ 30 km s ' of Model H 
are nearer the head. Model I displays a small range of velocity 
widths, which are also higher towards the head of the object. 
Some of the variation in the line widths is due to eddies being 
advected through the structure. 

5. DISCUSSION 

5.1. Champagne flow Model 

From Figure || we can derive approximate empirical expan- 
sion laws for the ionization front and shock into the neutral 
material for the photoevaporated champagne flow (Model A) 
and champagne flow plus stellar wind (Model C) cases. We 
find that for times 2xl0^<f<2xl0'' yrs a good approxi- 
mation to the ionization front velocity for both Models A and 
Cis 

y = 3.46f4'"kms-i, (3) 

while the shock front expansion for the same temporal range 
for Model A can be approximated by 

(4) 



v^A.15tf-^km^,-\ 



and for Model C by 



V = 4.99 f:" *"* km s" 



(5) 



where ti, is the time, f/10 yrs. If we work out the shock 
velocities for f = 2 x lO'* yrs, we find - 3.05 km s"' 




for Model A and = 3.20 km s ' for Model C. The shock 
velocity for Model C is slightly faster than that of Model A 
because the stellar wind shock has overtaken the ionization 
front shock. Although these results are specific to the stellar 
wind parameters and density distributio n (and time) , they are 
broadly consistent with observations by ^oshi et al (2005) of 
carbon recombination lines towards the ultracompact H II re- 
gion G35.20-1.74. These observations show velocities in the 
neutral gas in the PDR of 2.5 km s ' into the cloud, while the 
ionized gas has velocities away from the cloud, as expected 
for a photoevaporated flow (see ^enney et al. 2005 , for a de- 
tailed discussion of the velocities of the different components 
associated with ionization fronts). 

The ionization front velocity for both models at 20,000 yrs 
is 2.13 km s In these models, the front is almost D- 
critical and the ionized gas leaves the front at about half the 
sound speed {v ^ -5 km s"'), passing through a sonic point 
(v ^ -10 km s~ h by the time the ioni zation fraction becomes 

(gT ■ 



equal to unity (Penney et al. 2005). For more open (i.e 



less concave) ionization fronts the gas leaves the front even 
faster and accelerates more rapidly to the sound speed. Hence 



the net velocity of the ionized gas at the ionization front is 

3 km s"' in the rest frame of the cloud and becomes more 

blueshifted (in this scenario) as it flows away from the front. 
A popular misconception in the literature is that the velocities 
at the head of a champagne flow should be the same as those 
of the molecular cloud. 

We note here that the analytic expressions for the shape of 
the ioniza tion front i n a stro ng density gradient derived by, for 



example, [eke et al. (1980) are not correct since they do not 
take into account the motion of the gas, which has the effect 
of flattening the density gradient. Thus, the density gradient 
in the quasi-stationary champagne flow is not the same as the 
much steeper density gradient in the ambient material. The 
density gradient within these analytic ionized regions has no 
physical basis. In our numerical models, however, the hy- 
drodynamics and radiative transfer are solved simultaneously, 
giving a physical result. 



and collaborators (Fenorio-Tas 


rle 1979; iBodenheimer et al. 


1979; 


Bedijn & Tenorio-Taglel 


1981; Yorke et al. 1983), and 



more recently by Comeron (1997) are for the different physi 




cal setup of a spherical H II region breaking out from a dense, 
uniform cloud into a diffuse intercloud medium, where the 
distance of the star from the sharp cloud/intercloud interface 
is a determining parameter for the resultant flow. We feel that 
a more realistic density distribution would not be discontinu- 
ous, as in the above papers, but would show a smoother tran- 
sition between high cloud and low intercloud densities. Hy- 
drodynamical simulations of the formation and expansion of 
H II regions in disk-li ke clouds wit h strati fied density distri 



bution s were done by [Franco et al.| ( |1989| ) and [Franco et al 



(199C ), but in these models the star is placed at the midplane, 
which is a special place. For this reason, these models do not 
give rise to photoevaporated champagne flows. 

In the present work we adopted exponential density dis- 
tributions with different scale heights. These distributions 
have the advantage that the density depends only on the 
height, z- We also investigated a spherical density distribu- 
tion (Model F), where the density depends on both z and the 
cylindrical radius r. However, the results in this case are not 
qualitatively dissimilar from those of the exponential distri- 
bution with a large scale height (Model E). Models with large 
scale heights show the smallest departures from idealized H II 



regions, while small scale heights produce the greatest differ- 
ences. 

5.2. Bowshock Model 

In Figure ^ we compare our numerical results to the ana- 
lytical solution found by Wilkin (1996) for stellar wind bow- 
shocks in the thin-shell limit. The analytical solution for the 
shape of the bowshock is given by 



r{e) = ro CSC ^|2,{\ -6cot6) (6) 

where ro - (Mwi'm/^^Pai'iy^^ is the stand-off distance, M^,, v^, 
are the stellar wind mass loss rate and velocity respectively, 
is the ambient density and is the stellar velocity. The angle 
6 is measured with respect to the star from the symmetry axis. 
This model assumes that both the shocked stellar wind and 
the swept-up shell of ambient material are thin, i.e., that cool- 
ing is efficient in both shocked zones. The shape of the bow- 
shock should be compared to that of the contact discontinuity 
between the shocked stellar wind material and the swept-up 
ambient medium. 

The numerical result differs from the analytical result for a 
variety of reasons. The main reason is that the key assumption 
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Fig. 16. — Same as Fig. M but for Model D. 
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Fig. 18.^ Same as Fig. | but for Model F 
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Fig. 19.— Same as Fig. | but for Model G. 



E o 



HH-H- 



I I I I I 




I I I I I 



-10 10 




K i I 1 ^ 1 



I I I I I 



/A 




I I I I I 



I I I I I 



Offset Perpendicular to Axis 



I I I I I 



Hydrodynamics of cometary HII regions 



19 



E 



4-^ 



-H-+ 



10 



-10 



Offset Parallel to Axis 



Fig. 20.^ Same as Fig. | but for Model H. 
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Fig. 22. — Logaiithm of ionized number density for the simple bowshock 
model (Model B after 20,000 years nf nfnliitini]) with the analytic solution 
for the shape of the thin shell from Wilkin ( 199£ ) overplotted. 




2x10' 



3x10' 



4x10' 



Time [yrs] 



Fig. 23. — Distance of ionization front (solid line), stellar wind contact 
discontinuity (dashed line), and outer shock in neutral medium (dotted line) 
from the star as a function of time for the bowshock model (Model B). The 
dot-dashed horizontal line represents the position of ram-pressure balance 
between the stellar wind and the ambient medium at the apex of the bowshock 
formed by the moving star. 



of the analytical result is that cooling in the shocked stellar 
wind is instantaneous, making this region thin. In the present 
models the stellar wind shocks and a region of high temper- 
ature forms between the free-flowing wind and the swept-up 
ambient medium. The thickness of this subsonic region varies 
as the evolution proceeds since pressure balance across the 
contact discontinuity means that it is sensitive to what is hap- 
pening behind the outer shock. 

Additionally, in the numerical model the ionizing pho- 
tons and stellar wind switch on at the same time. Thus, 
a Stromgren sphere forms around the star almost instanta- 
neously, while the stellar wind bubble takes much longer to 



form. The stellar wind sweeps up a shell of material, which 
then cools on a timescale dependent on the postshock tem- 
perature and density. As the hydrodynamic flow develops, 
the initial photoionized region starts to expand, preceded by 
a shock wave, that sets the neutral gas ahead of the photoion- 
ized region in motion. It is into this moving medium that the 
stellar wind bubble expands. The impact of the hypersonic 
stellar wind on the surrounding medium set s up the familiar 
two shock pattern ( fPyson & Williams 1997), where the outer 
shock sweeps up the ambient medium. This swept-up shell 
forms initially within the photoionized region. However, due 
to the increased opacity, the ionization front retreats until it 
is contained with the shell. At this point, the shock ahead of 
the original ionization front is still expanding into the ambi- 
ent medium, although the previously ionized gas outside the 
shocked shell now recombines. The swept-up shell expands 
into this moving neutral gas and thickens. Instabilities are 
formed in the cooling shell, which lead to turbulent motions 
in the subsonic region between the outer shock and the in- 
ner stellar wind shock. The shell is initially thick, and only 
starts to become thinner once the outer stellar wind shock has 
overtaken the shock produced by the initial ionization front. 
The ionization front moves outwards within the stellar wind 
shell and the pressure of the photoionized gas is an additional 
confinement on the shocked stellar wind. Pressure gradients 
across the bowshock structure and transverse to it mean that 
the ram pressure balance across thin shells is inadequate to 
describe the evolution of the structure at these early times 
(t < 40, 000 vrs). 

In Figure 23 the distance of the ionization front and stellar 
wind contact discontinuity and the outer shock in the frame of 
the star are plotted for the bowshock model. We also indicate 
in this figure the position of theoretical ram pressure balance. 
In this figure we see how the initial ionization front retreats 
into the swept-up shell, while the shock that was sent out into 
the neutral gas travels for quite a time at roughly the stellar ve- 
locity (20 km s"') before being finally overtaken by the outer 
stellar wind shock. Once the ionization front is trapped in 
the shell, it remains there for the rest of the evolution. The 
noise in the position of the contact discontinuity occurs when 
instabilities form in the subsonic region between the shocks. 
The first peak in the positions of the ionization front and con- 
tact discontinuity, at about t ~ 4000 yrs, corresponds to the 
first peak in Figure EJ, interpreted as the first large eddy to be 
shed by the bowshock. The thickness of the shocked stellar 
wind region changes in response to pressure changes across 
the subsonic region, either due to the formation of instabili- 
ties or dynamics in the photoionized part of the shell. 

Interestingly, the position of the contact discontinuity in 
Figure ^ is settling down (after 40,0 00 yrs) to a value close 
to a! 1.5 ro, which was predicted by van Buren & McCray 



( 1988 ) from a simple model approximating the dynamics in 
the shocked wind by sonic flow in a hemispheric region be- 
tween the stellar wind shock an d the contact discontinuity. 
However, van Buren et al. (1990) argued that cooling via con- 
duction in the shocked stellar wmd would be so rapid that this 
layer would be very thin and its dynamics governed by mo- 
mentum conservation. We do not include thermal conduction 
in the present models but si ispect that it will not p lay as im- 
portant a role as invoked by van Buren et al. ( 1990 ). The im- 
portance of conduction fronts for the cooling of stellar wind 
bubbles has been cast into doubt by the discovery of diffuse 
X-ray emission from high-mass star-forming regions, which 
suggests that the shocked stellar winds inside H II regions 
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do not cool efficiently ( Townsley et al. 2005). It is probable 
that conduction is greatly reduced by the presence of mag- 
netic fields. Our results show that instabilities at early times 
in the shocked wind layer do enhance cooling and affect its 
thickness. However, this is only a transient feature, and once 
the initial instabilities have been advected away from the bow- 
shock head, the shocked stellar wind region is able to expand. 

5.3. Champagne flow plus stellar motion 

The density gradient plays a key role as to whether 
the dominant flow characteristic is champagne-flow-like or 
bowshock-like. In the steeper density gradient (Models G and 
H), a champagne flow sets up first, within which the stellar 
wind bubble develops, as for Model C. However, stellar mo- 
tion ensures that the outer stellar wind shock rapidly over- 
takes the shock ahead of the ionization front. The redshifted 
velocity of the ionized gas ahead of the star is less than it 
would be for a bowshock in a neutral medium because of the 
photoevaporated flow set up by the density gradient, which is 
blueshifted. The density gradient in the neutral material ahead 
of the moving star increases the ram pressure acting on the 
stellar wind and the photoionized shell. Both the shocked stel- 
lar wind zone and the ionized shell become thinner ahead of 
the star as it travels up the density gradient, with the shocked 
stellar wind zone almost disappearing. This is more pro- 
nounced in Model H, where the faster-moving star reaches 
higher density regions than the slower star of Model G. Be- 
cause the shocked stellar wind region ahead of the star be- 
comes extremely thin, it is essentially the ram pressure of the 
free-flowing wind which balances that of the stellar motion 
through the ambient medium. At the sides of the structure, 
a champagne flow is still in evidence and the pressure of the 
photoionized gas confines the stellar wind (see Figure |l2|). 

In the shallow density gradient (Model 1), on the other hand, 
the swept-up shell formed by the stellar wind confines the ion- 
ization front before a champagne flow has chance to form. 
The outer stellar wind shock soon overtakes the shock ahead 
of the ionization front. The stellar wind bubble structure ex- 
pands laterally downstream where the densities are lower and 
the swept-up shell becomes quite thick here. However, this 



never becomes a champagne flow (see. Figure 13 ). Ahead of 
the star the shocked stellar wind zone becomes thinner as the 
star advances up the density gradient, until it almost disap- 
pears and the ram pressure of the free-flowing wind balances 
that of the stellar motion through the ambient medium. The 
photoionized region ahead of the star is quite thick. 

Although the stellar velocity in Model G is only 5 km s \ 
the projected velocity range shown in Figure 19| is ~ 
15 km s with the velocities becoming redshifted towards 
the tail where the champagne flow is important. The actual 
velocities at the head of the object are slightly redshifted com- 
pared to the background cloud, but the projected velocities at 
the head of the object are blueshifted, unlike in the simple 
bowshock case (see Figure This shows that the pressure 
gradients within the shell due to the density distribution play 
an important role in advecting photoionized material in the 
shell at the head of the object round to the sides. Also, in 
projection we are seeing the near part of the champagne flow 
component, which has a strong influence on the projected ve- 
locities. The linewidths of this model also rise sharply to- 
wards the tail, reaching 40 km s"', suggestive of champagne 
flow. Model H also presents a projected velocity range of 
~ 15 km s"' for a stellar speed of 10 km s"'. The peak red- 
shifted velocities in this model are just downstream of the 



stellar position, and the peak linewidths are just ahead of it, 
which is more in keeping with a bowshock model. On the 
other hand, the velocities ahead of the star are blueshifted, 
as in Figure |9| for Model G, and this must be due to the in- 
fluence of the density gradient. The velocities then become 
blueshifted again far down the tail. 

The stellar speed in Model 1 is 10 km s ' and the ionized 
gas velocities in the model ahead of the star are low and red- 
shifted. However, even the shallow density gradient in this 
model ensures that the projected velocities ahead of the stel- 
lar position are blueshifted (see Figure |T]) and become red- 
shifted in the tail before turning blueshifted again. The full 
velocity range for this model is only about 8 km s"', princi- 
pally due to the shallow density gradient. The linewidths are 
highest around the projected stellar position. 

These combined models show that stellar motion in a den- 
sity gradient has important diff'erences to the standard bow- 
shock model. In particular, the velocity range is larger than 
that for the case of stellar motion only (e.g.. Model B) and the 
radial velocities change sign along the length of the object. 
These models are not hugely different from the champagne 
flow plus stellar wind cases, except that since there is effec- 
tively no shocked stellar wind zone ahead of the star (at least 
in the cases we have studied), then the peak emission occurs 
closer to star (see Figures landP). 

The intention of Models G, Hand I is to investigate what 
happens to the dynamics when density gradients, stellar winds 
and stellar motion are combined. In particular, we have not 
been overly concerned that the stars in our models start from 
low density regions. However, we could postulate a scenario 
in which stars, under the influence of the gravitational field 
of a molecular cloud, move in and out of local density peaks, 
and that our models represent a particular stage in this picture. 
Another possibility is that an external shock triggers star for- 
mation in the outer parts of a cloud and that, once formed, the 
massive star continues to coast in towards the cloud center, 
or is attracted to the center by the gravity of the main cloud 
core. In any case, the low density material is quickly advected 
off grid and replaced by either photoevaporated flow material 
(Models G and H) or material advected round the photoion- 
ized shell from the region of the head by the pressure gradient 
(Model I). In either case, the details of the initial environment 
are quickly lost by the model as it advances up the density 
gradient. 



Ffiinco et iil. ( p005[ ) have recently studied the evolution of 
H II regions in dense, structured molecular clouds where the 
star starts from a dense core and travels down the density gra- 
dient. They do not include the effect of stellar winds or discuss 
the dynamics of the ionized gas. The velocities of the stars 
vary between 2 and 12 km s"' (i.e., subsonic and supersonic 
with respect to the ionized gas). It is found that for slow mov- 
ing stars a champagne flow dominates, while for the faster 
moving stars a "bow" shape is formed when the star overtakes 
the ionized gas in expansion, since the gas can only expand at 
the sound speed and the star is moving supersonic ally. In the 
high density regions the star has left behind the gas recom- 
bines. This is an interesting scenario but unfortunately the 
authors do not include any analysis of the kinematics. It is 
probable that the inclusion of a stellar wind will modify the 
faster case, since the stellar velocity is not supersonic with 
respect to the stellar wind and the direction of motion is to- 
wards decreasing densities, so the ram pressure of the ambient 
medium is decreasing (or constant). In our models, the ram 
pressure increases as the star moves up the density gradient 
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Fig. 24. — Time evolution of projected velocities at half peak flux (solid 
line), quarter peak flux (dashed line) and one eighth peak flux (dotted line) 
for inclination angle = 45°: Top panel — Model B (bowshock); Bottom 
panel — Model C (champagne flow plus stellar wind). 



and so it is quite possible that the star overtakes the swept-up 
wind shell with its trapped ionization front if the evolution is 
followed for long enough. 

5.4. Effect of Eddies 

In Figure ^ we examine the evolution of the line-of-sight 
velocity at half, quarter and one eighth of the peak flux along 
the symmetry axis of the numerical model for an inclination 
angle of / = 45°. The large peaks in these curves correspond 
to eddies being shed from the head of the cometary H II re- 
gion and travelling downstream, hence there is a phase lag 
between peaks on each line. In the champagne flow plus stel- 
lar wind case an "average" profile can be discerned for each 
curve. However, in the bowshock case the curves are far too 
noisy. It is not possible to draw any conclusions from these 



plots to aid as a diagnostic tool for discriminating between 
bowshock and champagne flow plus stellar wind models. We 
note that in the bowshock case there appears to be a slight 
trend for the velocities further down the slit to rise with time, 
while in the champagne flow plus stellar wind case there is 
a trend for all velocities to become slowly more blueshifted 
with time. 

The spectra obtained by, e.g., f^^umsden & Hoare ( 1996 



1999) are not as noisy as our theoretical spectra. This could 



partially be due to the resolution of the observations, but it 
also suggests that the eddies we see forming in our simula- 
tions are dissipated more quickly in reality. In order to in- 
vestigate this possibility numerically, we would have to in- 
corpor ate a turbulence model, such as that outlined by Falle 
( |1994| ), in our numerical scheme. This is beyond the scope or 
the present paper. 

5.5. Comparison with observations 
5.5.1. G29.96-0.02 

G29.92-0.02 is, perhaps, the most well studied of all 
cometary ultracompa c t H II regions. It was first observed by 
Wood & Churchwell (1989) at radio wavelengths. The mor- 
phology of this region is that of a bright arc of emission at 
the head trailing off into a low-surface brightness tail. The 
background cloud systemic velocity is best determined from 
molecular line ob servations, which trace t he bulk of the gas . 
Both single dish ( phurchwell et al.l |1992|; [kim & Kool TO03| ) 
and interferometric ( ^'atap et 

servations of optically thin CO and CS lines give velocities of 
97-98 km s-i. 

The infrared observations of Lumsden & Hoarq (|l996l 
1999| ), [Martm-Hernandez et al.| ( |2003D and |Zhu et al "| (POOSD 



provide the most detailed kinematic results of G29.96-0.02 
with which to compare our numerical simulations. Lu ms- 
den & Hoare'- ^TWf] ) and [Martfn-Hernandez et al.| ( |2(j53| ) ob- 
tained long slit spectra of t he Bry hydrogen recombination 



line, while Zhu et al. (2005) obtained mid-infrared observa- 
tions of the [Ne II] 12.8 fim fine-structure line. These fine- 
structure line studies are particularly promising because the 
FWHM line width due to thermal broadening of Ne^ at 10"^ K 
is only ~ 4.8 km s"', whereas that of H^ is ~ 21 km s"\ 
hence the kinematics of the ionized gas can be separated from 
the thermal and turbulent broadening contributions if the ve- 
locity resolution of the observations is good enough (as long 
as the Ne^ ion is present in the whole volume). 

These detailed observational studies find a radial velocity 
range of 15 to — hlO km s"' for the ionized gas with re- 
spect to the molecular cloud (assuming a background cloud 
velocity of 98 km s"^), with the highest blue-shifted veloci- 
ties ahead of the emission peak of the object, red-shifted ve- 
locities behind the emission peak and blue-shifted velocities 
again in the tail. The line widths of the Bry data are highest 
at the head, as are the [Ne II] line widths. In the past, this 
object has been modeled as a bowshock H II region, with an 
inferred stellar motion of 20 km s~' in order to give the ve - 



locity range (e.g.. Afflerbach et al 



1996 



Martm-Hernandez et al 



1994: Lumsden & Hoare 



2002; Zhu et al. 2005). How 



ever, for a slit placed along the symmetry axis of the cometary 
H II region, simple bowshock models give either wholly red- 
shifted or wholly blueshifted velocities with respect to the 
background cloud, depending on the angle of inclination to 
the line of sight. For this reason, several authors (L umsden 
& Hoare |1996^ |1999| ; [Martm-Hernandez et al.|p003| ; Zhu et 
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al. 2005 ) have argued that a simple bowshock model cannot 
reproduce the observed velocity behavior. Although a stellar 
velocity of 20 km s"' can give the observed velocity range, it 
cannot reproduce the velocity offset with respect to the molec- 
ular cloud, which requires a lower stellar velocity. On the 
other hand, simple champagne flow models would not give a 
limb-brightened mo rphology, nor reproduce t he velocity and 
linewidth behavior. Lumsden & Hoare ( 1996 ) suggested that 
stellar winds and density gradie nts might be impo rtant for the 



kinematics of this object, while Zhu et al, (12005) argued that 
pressure gradients within the parabolic shell are needed to in- 
crease the velocity range. 

Our numerical models comprising a star with a stellar wind 
in a density gradient (e.g.. Models C, G, H and I) can all re- 
produce the blueshifted to redshifted to blueshifted velocity 
behavior exhibited by G29.96-0.02, together with a strong, 
arc-like emission at the head, assuming an inclination angle of 
45° to the line of sight. However, the maximum velocity range 
our models achieve is only ~ 15 km s ' before the flux drops 
below 1% of the peak value, as against the ~ 20 to 25 km s"' 
seen in the observations. The velocity range increases as the 
density gradient becomes steeper, so it is conceivable that a 
steeper gradient than the ones we have investigated would 
give the observed velocity range. The spatial variation of the 
flux of Bry al ong the symmetry axis in the hig h resolution ob- 
servations of Martm-Hernandez et al. ( 2003 ) bears a striking 
resemblance to that of our Model I, as does the variation of the 
velocity centroid of the same line along the slit (ignoring the 
blue-shifted shoulder seen in the observations, which appears 
to be unconnected with the H II region and is probably due 
to motions associated with the nearby hot core source). Even 
the "knee" in th e observed jBry fl ux, also present in the [Ne II] 
observations of Zhu et al. (2005), is seen in the model (where 
it is due to the line of sight cutting through thick, dense re- 
gions of ionized shell on both sides of the star), although this 
is serendipitous. Models C, G and H, however, have a larger 
velocity range than Model I. The linewidths in Models H and 
I are highest at the head of the object for this slit, although 
not ahead of the pe a k emission, as se en in the observations by 
|l.umsden & Hoare| (1996 , 1999) and Martm-Hernandez et al. 
(2003). The large linewidths seen in Model 1 are most likely 



due to instabilities in the photoionized shell and are thus tran- 
sient. The linewidths in Model H, on the other hand, appear 
to be the result of the general flow pattern. 

The H2 observed in this object comes from the region out- 
side of the ionization front. However, it is not likely to come 
from the quiescent background cloud, rather it arises in the 
region between the shock ahead of the ionization front in the 
neutral medium and the ionization front, probably due to flu- 
orescence, and thus will have a red-shifted velocity with re- 
spect to the cloud and the io nized gas. This is suggested by 
the observations reported by Martm-Hernandez et al. (2003), 
where the H2 is redshifted by ~ 2 km s"' with respect to the 
background cloud velocity as determined from molecular line 
observations. Our models show that a shock front advances 
slowly into the neutral medium with velocities of order a few 
km s"', which is consistent with this velocity redshift. 

The parameters of the numerical models have not been tai- 
lored to fit G29.96-0.02, but the qualitative agreement is en- 
couraging. We have shown that models including a strong 
stellar wind and a density gradient are able to reproduce the 
morphology, velocity behavior, and to some extent the line- 
width behavior of this object. Stellar motion appears to mod- 
ify the line-width behavior of an essentially champagne-type 



flow and so we suggest that G29.96-0.02 could be the H II 
region formed by a star with a strong stellar wind in a density 
gradient, possibly with the star moving up the gradient. 

5.5.2. Hydrogen radio recombination line studies 

A detailed examination of H92q' radio recombination line 
ob servations of three cometary H II regions was presented 
by |Garay et"ar (1994). They found velocity ranges of up to 
12 km s"' between the head and tail regions of these objects. 
Two of these objects were interpreted as champagne flows, 
mainly because the velocities at the head were approximately 
equal to the background cloud velocities as measured from 
observations of NH3. The third object was interpreted as a 
bowshock flow since the velocities were all redshifted with 
respect to the background cloud. However, even in a sim- 
ple photoevaporated champagne flow model such as that pre- 
sented in Model A, we would expect to find a velocity shift 
between the ionized gas at peak emission and the background 
cloud (unless the object is in the plane of the sky), since the 
gas leaves the ionization front with a velocity about half the 
sound speed and then accelerates (see, e.g.. Figure ||). Also, in 
the champagne flow objects, the linewidth variations between 
head and tail are not entirely consistent with the champagne 
flow interpretation: in one object the linewidths decrease to- 
wards the tail and in t he other they are roughly constant. 



Uaray et al. (1994) discuss the possibility that the stellar 
wind modifies the flow, but their analytic treatment requires 
the star to be at the center of a radially decreasing density 
distribution. In our models we are able to place the star and 
its stellar wind at arbitrary points in density distributions with 
steep and shallow gradients. We find that the combination of 
a stellar wind and a density gradient changes the linewidth 
behavior, since the ionized flow is forced to flow around the 
outside of the stellar wind b ubble (e.g.. Models C, E and F). 

Cyganowski et ST ( 2003 ) have recently reported hydrogen 
recombination line radio observations of dual cometary H II 
regions in DR 21. The radial velocities of the ionized gas 
suggest a bowshock interpretation since they are higher at the 
head. However, for one of the H II regions, it is reported else- 
where ( Roelfsema et al. 1989 ) that velocities down the tail be- 
come increasingly redshifted with respect to the background 
cloud, while velocities at the head are blueshifted. Further- 
more, the linewidths for both objects suggest that there is an 
increase towards the tail, which is not expected for a bow- 
shock model. Models incorporating stellar winds and density 
gradients, such as our Models C, F and G, can demonstrate 
this spatial variation in velocity and linewidth . The "hybrid" 
picture developed by Cyganowski et al. ( 2003 ) to explain the 
motions they see in DR 21 has some similarities with our 
models. 

From observations of the neutral, atomic and molecular 
;as associated with t he cometary H II region Gl 1 1.67+0.37, 
Lebron et al, (2001) find that both the neutral and ionized 



gas exhibit velocity gradients from the head to the tail of the 
cometary region, but that the neutral gas has a lower velocity 
range than the ionized gas and is detected principally towards 
the tail where the PDR is thickest. Our numerical models also 
demonstrate this behavior Although we do not show results 
for the neutral gas here, it will be the subject of a future paper 

6. CONCLUSIONS 

We have performed numerical simulations of cometary H II 
regions forming in media with density gradients and where 
the ionizing star may additionally have either or both a strong 
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stellar wind and a supersonic motion with respect to the neu- 
tral gas. Analysis of our results and a comparison with obser- 
vations from the literature show that: 

1. Photoevaporated champagne flow models with stellar 
winds can have limb-brightened morphologies, which 
simple champagne flows cannot. This is because the 
stellar wind bubble acts as an obstacle to the photoe- 
vaporated flow, which must divert around it. 

2. Depending on the steepness of the density gradient, the 
morphologies of the champagne flow plus stellar wind 
models can appear more limb-brightened (steep gradi- 
ent) or more sheU-Uke (shaUow density gradient). 

3. In the absence of thermal conduction, numerical bow- 
shock simulations (in uniform media) are not well de- 
scribed by momentum-conserving, ram pressure bal- 
ance analytical models. This is because the shocked 
stellar wind does not cool efficiently and its thermal 
pressure must also be taken into account. 

4. In both champagne type models and bowshock models 
a shock runs ahead of the ionization front in the neutral 
material. Typical velocities in the neutral gas behind 
these shocks are a few km s"^ These velocities could 
influence H2 and carbon recombination line measure- 
ments of the systemic cloud velocity, since this emis- 
sion is most likely to come from the photon dominated 
region, which is a thin region just outside the ionization 
front and thus whoUy contained in the zone behind the 
shock in the neutral gas. 

5. Models in which a star with a strong stellar wind 
advances supersonically up a density gradient are 
champagne-flow dominated if the gradient is steep 
enough and stellar wind shell dominated if the den- 
sity gradient is shaUow. These models produce Umb- 
brightened H II regions. In these models, the shocked 
stellar wind region is extremely thin due to the increas- 
ing ram pressure as the star advances up the density 
gradient, and the fact that the shocked wind material 
is quickly advected away from the region of the head 
due to the density gradient. 

6. Models of stars with stellar winds moving in density 
gradients result in the peak emission being closer to the 
star than for similar models without stellar motion be- 
cause the shocked stellar wind zone is extremely thin 
and so the photoionized shell is closer to the star. 

7. The radial velocity range of ~ 20 km s ' seen in obser- 
vations, coupled with a limb-brightened morphology. 



Afflerbach, A., Churchwell, E., Hofner, P., & Kurtz, S. 1994, ApJ, 437, 697 

Arquilla, R., & Goldsmith, P. F. 1985, ApJ, 297, 436 

Arthur, S. J., & Falle, S. A. E. G. 1993, MNRAS, 261, 681 

Arthur, S. J., Kurtz, S. E., Franco, J., & Albarran, M. Y. 2004, ApJ, 608, 282 

Bedijn, P. J., & Tenorio-Tagle, G. 1981, A&A, 98, 85 

Bodenheimer, P, Tenorio-Tagle, G., & Yorke, H. W. 1979, ApJ, 233, 85 

Chevalier, R. A., & Clegg, A. W. 1985, Nature, 317, 44 

Churchwell, E., Walmsley, C. M., & Wood, D. O. S. 1992, A&A, 253, 541 

Comeron, F 1997, A&A, 326, 1 195 

Comeron, F, & Kaper, L. 1998, A&A, 338, 273 

Cyganowski, C. J., Reid, M. J., Fish, V. L., & Ho, P T. R 2003, ApJ, 596, 344 
De Pree, C. G., Wilner, D. J., Deblasio, J., Mercer, A. J., & Davis, L. E. 2005, 
ApJ, 624, LlOl 



can be reproduced either by a fast-moving ~ 20 km s" 
star or by a star with a stellar wind in a steep density 
gradient. The key observational discriminant is the ra- 
dial velocity pattern along the object with respect to the 
ambient material. For an inchnation angle of i = 45° 
to the line of sight, for example, the velocities of the 
bowshock model will all be redshifted with respect to 
the ambient cloud, whereas the champagne flow plus 
stellar wind model gives blueshifted velocities at the 
head, redshifted velocities around the stellar position, 
and blueshifted velocities again in the tail. Similar ve- 
locity behavior is seen in the models that also include 
stellar motion, showing that for the modest stellar ve- 
locities considered (5 and 10 km s '), the density gra- 
dient is more important than the star's speed for deter- 
mining the velocities in the photoionized shell. 

8. The linewidths in simple champagne flow models are 
highest in the tail of the cometary H II region, while 
those of bowshock models are highest ahead of the 
star. In our models combining density gradients, stellar 
winds and stellar motion we find that the champagne 
flow dominated flows (slow-moving star in steep den- 
sity gradient) have larger velocity widths towards the 
tail, while for faster moving stars in shallow density 
gradients the largest line widths occur further up the 
cometary structure at the stellar position. Such inter- 
mediate hnewidth information is also reported in ob- 
servations. 

Our results show that the most likely explanation for the ob- 
served morphologies and kinematics of cometary compact 
H II regions is the outcome of ionizing stars with strong stel- 
lar winds in density gradients, possibly with the star moving 
up the gradient. It is important to obtain a reliable determi- 
nation of the systemic cloud velocity in order to be able to 
discriminate between possible models. 
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